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ABSTRACT 

We define and analyze the photometric orbit (PhO) of an extrasolar planet 
observed in reflected light. In our definition, the PhO is a Keplerian entity with 
six parameters: semimajor axis, eccentricity, mean anomaly at some particular 
time, argument of periastron, inclination angle, and effective radius, which is 
the square root of the geometric albedo times the planetary radius. Prelimi- 
narily, we assume a Lambertian phase function. We study in detail the case of 
short-period giant planets (SPGPs) and observational parameters relevant to the 
Kepler mission: 20 ppm photometry with normal errors, 6.5 hour cadence, and 
three-year duration. We define a relevant "planetary population of interest" in 
terms of probability distributions of the PhO parameters. We perform Monte 
Carlo experiments to estimate the ability to detect planets and to recover PhO 
parameters from light curves. We calibrate the completeness of a periodogram 
search technique, and find structure caused by degeneracy. We recover full or- 
bital solutions from synthetic Kepler data sets and estimate the median errors in 
recovered PhO parameters. We treat in depth a case of a Jupiter body-double. 
For the stated assumptions, we find that Kepler should obtain orbital solutions 
for many of the 100-760 SPGP that Jenkins & Doyle (2003) estimate Kepler will 
discover. Because most or all of these discoveries will be followed up by ground- 
based radial-velocity observations, the estimates of inclination angle from the 
PhO may enable the calculation of true companion masses: Kepler photometry 
may break the "msini" degeneracy. PhO observations may be difficult. There is 
uncertainty about how low the albedos of SPGPs actually are, about their phase 
functions, and about a possible noise floor due to systematic errors from instru- 
mental and stellar sources. Nevertheless, simple detection of SPGPs in reflected 
light should be robust in the regime of Kepler photometry, and estimates of all 
six orbital parameters may be feasible in at least a subset of cases. 



Subject headings: techniques: radial velocities — techniques: photometric — plane- 
tary systems 
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1. Introduction 

We address Kepler's observations of "photometric orbits" (PhOs), which is our term for 
the variable flux of reflected starlight from short-period giant planets (SPGPs). We believe 
that an improved understanding of this enigmatic class of objects will be an important legacy 
of the Kepler mission. We have conducted a ministudy of Kepler's expected performance 
in detecting SPGPs and in estimating their orbital parameters via the PhO. We have found 
unexpected scientific benefits. 

A "Keplerian data set" is an observational record governed by the classical motion of 
two bodies interacting by the force of gravity. Radial-velocity, astrometric, and photometric 
observations of extrasolar planets are prime examples. In fact, the three types provide 
overlapping and complementary information, which creates synergism when more than one 
data type is available, as we can expect in many cases for SPGPs discovered by Kepler — at 
least for radial velocity and photometry (including both transits and reflected light). 

When first discovered in the 1990s, SPGPs were a surprise, and we still do not under- 
stand their basic properties. Why do they exist in such abundance? How were they created? 
After creation, did they migrate in distance from the star, and if so, by what mechanism? 
Why do the orbits of many SPGPs have significant eccentricity, if tides should efficiently 
circularize them? (See Fig. 1.) What are the optical radii and actual masses of SPGPs, 
and what are their physical conditions? Of what types of stellar and planetary systems are 
SPGPs members? 

In the next decade, new information from a variety of sources will help address these 
questions, and we expect that Kepler's observations will be important. Kepler's discoveries 
will increase the number of known SPGPs, and its measurements of PhOs will increase the 
value of other data sets by independent estimates of some planetary parameters, and will 
provide some information not otherwise available. 

The current Kepler science team plans to search for SPGPs by seeking peaks in the 
periodograms (power spectra) of stellar light curves and to determine geometric albedos 
for transiting SPGPs (Jenkins & Doyle 2004). This paper goes well beyond those modest 
ambitions. 

We analyzed many simulated Kepler data sets to explore how reflected-light observations 
might provide additional information about SPGPs. (Fig. 2 shows one example, treated 
in depth in §8.) These Monte Carlo experiments allow us (1) to provide a preliminary 
calibration of the completeness of the periodogram search technique. Such calibration is 
generally informative and will enable estimates of the true occurrence frequency of SPGPs 
as a function of planetary parameters. A sufficiently strong signal allows us (2) to provide 
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complete solutions of the PhO. So far as we can detect, no one else has previously thought 
this possible. 

The benefits of full PhO solutions include: (1) using the estimated orbital inclina- 
tion to determine the true planetary mass when radial-velocity measurements are available; 
(2) improved albedo estimates for eccentric transiting planets, by allowing data points to 
be interpreted with their correct planetary phase angles; (3) independent estimates of the 
four orbital parameters in common with radial velocity; and (4) possible estimates of the 
planetary phase function in very high signal-to-noise cases. 

PhO observations may be difficult. There is uncertainty about how low the albedos of 
SPGPs actually are, about their phase functions, and about a possible noise floor due to 
systematic errors from instrumental and stellar sources. Nevertheless, simple detection of 
SPGPs in reflected light should be robust in the regime of Kepler photometry, and estimates 
of all six orbital parameters may be feasible in at least a subset of cases. 



2. Theory of the PhO 

To the best of our knowledge, a theory of the PhO does not exist in the literature, and 
so we provide a summary account. 

The most convenient observable of the PhO is based on — but not exactly equal to, as 
we will see — the ratio of planetary to stellar flux, which is the function: 

h=^m , (i) 

where r is the radial distance between the star and planet; R e s is the effective planetary 
radius 

R cS = ^R p , (2) 

where p is the geometric albedo, and R p is the planetary radius; and is the phase 
function, where (3 is the planetary phase angle, which is the planetocentric angle between 
the star and observer. Nominal values of p are 0.66 for a conservative Lambert sphere (a 
ping-pong ball) and 0.5 for Jupiter at visible wavelengths. Thus, for Jupiter, R cS = 7.9 R®, 
in terms of the Earth radius. 

We expect to decrease monotonically from 1 to over the range of /3 from to 
7r. In the absence of better information, we adopt the Lambertian phase function in the 
ministudy: 

= sin/? + (7r-/?) cos/3 

7T 
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Fig. 1. — Eccentricities and periods of the 40 extrasolar planets with period T < 10 
days, as listed in the catalogue at the California & Carnegie Planet Search webpage 
(http://exoplanets.org/) on March 23, 2006. 



Table 1. The ministudy's "planets of interest" 



Parameter 


Minimum 


Maximum 




Distribution 


Semimajor axis (a) 


0.03 


0.09 


AU 


uniform 


Eccentricity (e) 


0.0 


0.5 




uniform 


Initial mean anomaly (Mo) 





2tt 


radians 


uniform 


Argument of periastron (u) 





2tt 


radians 


uniform 


Inclination angle (i) 





7T 

2 


radians 


uniform on sphere 


Effective planetary radius (R c s) 





20 




uniform 


Period (T) 


1.9 


9.9 


days 


computed from a 



Note. - The random deviate for i is cos 1 (1-72.), where TZ produces pseudorandom 
numbers in the range 0-1. 
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Fig. 2. — Simulated data set created for the Jupiter body-twin with 20 ppm noise. The mean 
has been subtracted, and the data is folded on the orbital period. The curves show the true 
(red) and best-estimate (blue) solutions. We treat this data set in depth in §8. Although 
one's first impression might be that this is a very noisy data set, in fact the signal-to-noise 
ratio is S/N = 48. (S/N is discussed in §10.) The planet is easily detected by Kepler with 
the assumed performance, and all orbital parameters are robustly estimated (although there 
are systematic effects due to uncertainty in the phase function; see §9). 
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With superior data, the opportunity may arise to introduce into the analysis, and to estimate, 
one or more additional PhO parameters associated with the phase function. In §10 we discuss 
the systematic errors introduced when a data set is analyzed with the wrong phase function. 

As usual, the planetary position relative to the star is fully described by seven parame- 
ters: a (0 — oo), the semimajor axis of the orbit with a focus at the star; e (0 — 1), the orbital 
eccentricity; M (0 — 2rr), the mean anomaly at an arbitrary ti, taken here to be the time of 
the first observation; u (0 — 2ir), the argument of periastron; % (0 — 7r), the inclination angle; 
Q (0 — 2tt), the position angle of the ascending node; and T (0 — oo), the orbital period. 

Kepler's photometric observations are spatially unresolved and therefore do not con- 
strain Q, which we set to zero without loss of generality. Similarly, i is nondegenerate only 
in the range — ir/2. 

If we can estimate the sum of the planetary and stellar masses, m p and m s , then Kepler's 
Third Law, informed by Newton's Law of Gravity, provides a relationship between a and T: 

3 _ G(m p +m s ) T2 
a ~ 47r2 > W 

where G is the gravitational constant. We take advantage of Eq. (4) to retire T as an 
independent parameter by assuming that we can neglect m p relative to m s , and that we will 
have an estimate of m s available from the spectrophotometric characteristics of the star. 

To compute the planetary position in space, which is needed to compute r and (3 in 
Eq. (1), we establish a Cartesian coordinate system {x, y, z} with its origin on the star. 
The plane of the sky is the domain {x, y, 0}, with +x north, +y east, and +z in the 
direction of the observer. The base {a, e} orbit, with uo = and % — 0, lies in the plane of 
the sky, with periastron located on the +x axis, and the motion of the planet at periastron 
is in the +y or east direction. The base orbit is direct, meaning that position angle always 
increases with time. The first Eulerian angle, iv, is a right-hand rotation around +z axis, 
which moves periastron to the east. The second Eulerian angle, i, is a left-hand rotation 
around the +x axis. 

The planetary position in the base orbit at time t is 

{x b , y h , z h } = {r cosu, r sinv, 0} , (5) 

where the radial distance between the star and planet is 

„(!-«») 

(1 + ecosi/) ' W 
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where v is the true anomaly, which is the root of the equation 



v /(1 + e) E 
tan- = y ir — y tan- , (7) 

where E is the eccentric anomaly, the root of Kepler's Equation: 

E-esin E = M , (8) 

and M is the mean anomaly: 

M = M + 2n t - 1 ^ , (9) 
where M is the initial mean anomaly, and t\ is the time of the first observation. 
The planet's position in space is then 

{x, y, z) = r{ — sin v sin uj + cos v coscj, cos v cosi sin uj + sin v cosi cos uj, 

— cosu sin 2 sinz/ — cosu sini sincj} , (10) 



from which we can compute: 



/3 = arctan — — . (11) 

z 



The results r and (3 from Eqs. (6) and (11) permit the computation of h[t] in Eq. (1). 
§8 provides a particular, in-depth example of a PhO. 



3. Keplerian Data Sets of the PhO 

In the ministudy, we constructed each simulated data set from particular values for ten 
parameters, which served as the controls of our Monte Carlo experimentation. The parame- 
ters, grouped in three subsets, were generated as follows: (1) PhO parameters {a, e, M , iv, i, 
R c s} were drawn from appropriate random deviates over ranges relevant to Kepler's SPGPs 
(see Table 1). These ranges and deviates define "the planets of interest." (2) We chose 
the three observational parameters, {a, duration, cadence} = {2 x 10~ 5 , 3 years, 6.5 hours}, 
where we used the values of photometric uncertainty, a, and cadence provided in NASA Re- 
search Announcement NNH07ZDA001N-KPS, soliciting Kepler Mission Participating Sci- 
entists. This document states: "The mission will achieve a photometric precision of 20 
parts-per-million (ppm) on a m v = 12 magnitude G2V star, for a 6.5-hour integration and 
using differential ensemble processing." Our choice of the duration of the observations is 
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arbitrarily 75% of the planned mission duration. (3) With no particular star in mind, we 
chose m s — 1 M , in terms of the solar mass. 

Each simulated data set had the form 

{t, h} = {{tx, {t 2 , /i 2 },-{ijv, /ijv}} , (12) 

containing a number of data points iV equal to the integer part of ^^ a ^ Qn , or N = 4046. 

The values of fa were generated by a normal random deviate with mean value equal to the 
value of h in Eq. (1) at the value of t and with standard deviation a. In other words, the 
value of h was not smoothed over a 6.5 hour interval. 



4. Periodogram Planet Search 

At least for non-transiting SPGPs, Kepler's first indication of a planet will come from 
a periodogram analysis. Also, the information in the periodogram about the planetary 
period (semimajor axis) will provide the starting point for finding a full solution for the 
PhO. For these reasons, understanding the strengths and limitations of the periodogram 
is a cornerstone of PhO research. In addition, we recognize an opportunity to calibrate 
the periodogram in terms of search completeness, which will enable the estimation of the 
underlying SPGP population. 

The periodogram is the Fourier transform of the autocorrelation function. In our nor- 
malization, the estimated periodogram pg is the vector 

1 - 2 

Pglfreqj] = — G (freq 3 ) , (13) 
where G is the discrete Fourier transform of the data record: 

N 

/A" 

where 



1 N 

G\fre qj ] = ^=J2 h ^ ifreq > tk » ( 14 ) 



k=i 



* N cadence ' 

and j = 0, l,2...y indexes the non-redundant frequencies. The highest value, freqjy/2 = 

4 , is the Nyquist frequency of critical sampling. 

2 cadence J ^ ^ J ^ ° 

The ministudy's protocol was, first, to smooth the periodogram with a Daniell window 
of second order; second, to locate the highest point in the periodogram for j > and record 
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its frequency (pgFreq) and peak value (pgPeak); and third, to determine if the planet was 
"detected" by this test: 

rn 7 1 2 3 
pgFreq = -, -, or - , (16) 

within 2%. (For higher e, dominant spectral power may be transferred from the fundamental 
frequency to harmonics.) 

Figure 3 shows the parameter-specific photometric completeness (PSPhC) of the peri- 
odogram search technique. Here, we created 15 frames of 100 x 100 cells (1% resolution) to 
compile and display the results of the Monte Carlo experiment. This "paramagram" is an 
exhaustive two-dimensional projection of parameter space, the frames comprising all possi- 
ble pairs of PhO parameters. After determining which cell in each frame corresponded to 
the data set's control values of {a, e, M , cu, i, R e s}, we incremented the "number tested" 
(nTested) in that cell by one, and if the planet was detected, we also incremented the "num- 
ber detected" (nDetected) in that cell by one. After accumulating all results, we estimated 
PSPhC in each cell using the general definition: 

nDetected 

estimated completeness = — — — , (17) 

nTested 

which is the best estimate of the probability of detecting a planet with parameters in a cell's 
range (Table 1), assuming its presence. 

PSPhC is independent of M , because we are observing a large number of cycles. There- 
fore, any direct variation — or lack of variation — of PSPhC with respect to another parameter 
is seen best in the frame of the paramagram it shares with M . For example, frame 6 shows 
that the periodogram technique is neutral with respect to e over the range of e in Table 1. 

The planetary signal vanishes for face-on, circular orbits {e = 0, % = 0}, which is 
indicated by the point of low PSPhC in the bottom-left corner of frame 8. The streak of 
degraded PSPhC extending diagonally up and right from this corner is due to the ability 
of increasing i, when u m to cancel increasing signal near periastron from the reduced 
periastron distance a(l — e) due to e increasing. (Increasing i tilts the dark side of the planet 
toward the observer most efficiently at uo xs as documented by the horizontal streaks of 
degraded PSPhC in frames 3, 7, and 10, and the vertical streaks in frames 13, 14, and 15.) 
PSPhC is reduced because the signals of some of the planets touched by the effect, which 
were otherwise borderline-detectable, now fail to exceed the detection threshold. 

We expect the peak-to-trough amplitude of the PhO, peakToTrough = h max — h min , to 
correlate with (1) detectability via periodogram and (2) confidence in an estimated orbital so- 
lution. In other words, peakToTrough — a transcendental function of {a, e, M , cu, i, R e s} — is 
also a useful control parameter for our Monte Carlo experiments. That is, peakToTrough can 
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Fig. 3. — Ministudy results for the parameter- specific photometric completeness (PSPhC) of 
periodogram searches for extrasolar planets. As shown by the red lines in frame 5, PSPhC = 
10, 50, 90, and 99% is achieved for R eS /a « 35, 53, 91 & 161 R&/AXJ, respectively. This 
functional dependence is expected from Eqs. (1) and (6), which indicate the PhO is invariant 
under changes in R e s and a if the value of R e s/a is preserved. The magenta "J" indicates 
the Jupiter body- twin studied in §8. 
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be calculated when a synthetic data set is created and carried along to help interpret the 
results. From Eq. (1), the maximum value is peakToTrough = 8 x 10~ 4 for e = for the 
planetary population defined by Table 1. 

Figure 4 plots the 10 7 periodogram results versus peakToTrough and pgPeak. In green, 
we show the log probability density of detected planets parsed on a 1000 x 1000 grid. In 
the sinusoidal approximation of h-h min , the second-order power law evidenced in Figure 4, 
pgPeak oc peakToTrough 2 , is expected from Parseval's Theorem: 

N/2 1 N h I 2 N 

P9 [f rea j] = 4^ \ h k \ ~ 32^: PeakToTrough 2 , (18) 

j=0 k=l 

because 

N/2 

^2 pg [freqj] <x pgPeak , (19) 
j=o 

where the constant of proportionality depends only on the chosen smoothing of the peri- 
odogram. The minimum detectable value of peakToTrough is seen to be approximately 2a. 

Along a vertical at any value of peakToTrough in Figure 4, the cases with higher eccen- 
tricities are smeared downward from the upper-left edge of the probability density distribu- 
tion, due to increasing leakage of power from the fundamental frequency into harmonics as 
e increases. 

For Figure 4, search completeness was again computed via Eq. (17), but this time 
specifying peakToTrough and pgPeak, without regard to the PhO parameters. We call these 
quantities the amplitude-threshold photometric completeness (AThPhC) and peak-threshold 
photometric completeness (PThPhC), shown by the color-coded bands along the abscissa 
and ordinate. For example, we find AThPhC = 0.5 for peakToTrough = 4.5 x 10~ 6 and 
PThPhC = 0.5 for pgPeak = 2.5 x 10~ 10 . (The second number depends on the choice of 
smoothing, but the first number is universal, corresponding to S/N = 7; see discussion of 
S/N is §10.) 

The fraction of all 10 7 planets in the ministudy that were detected according to Eq. (16) 
is 0.91. This quantity is the estimated ensemble photometric completeness (EPhC), where 
the "ensemble" is the planetary population of interest defined in Table I. 

The benefits of studying the search completeness of extrasolar planetary observations 
are fourfold, at least (Brown 2004a, 2005, 2009). First, such study establishes realistic ex- 
pectations by clarifying search power in objective terms. Second, it provides a scientific 
metric to inform trades during mission development, and to optimize variable aspects of the 
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Fig. 4. — In green, the distribution of log-relative probability density for detected planets on 
the plane spanned by peakToTrough ("amplitude") and pgPeak ("peak"). The color-coded 
bands along the ordinate and abscissa show the completenesses PThPhC and AThPhC. The 
grayscale stripe illustrates periodogram-power leakage due to eccentricity. The red ocher line 
indicates the median value of pgPeak for detected planets, 2 x 10~ 8 . In §7, we discuss the 
universe of full orbital solutions for the population of detected planets above this line. In 
§8, we discuss the details of one particular planet, the Jupiter body-twin indicated by the 
magenta "J". 
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observations in the operational period. Third, completeness permits the underlying plane- 
tary population to be estimated, by taking into account intrinsically equivalent planets that 
happen to be undetectable because of unfavorable extrinsic characteristics. (We learn where 
planets can hide, how many are hiding, and compensate statistically.) Fourth, complete- 
ness studies offer unique insights into the complexities of Keplerian observations, as well as 
glimpses of their intellectual beauty — and occasional obscurities. 



5. The Observable h-<h> 

Rather than h, Kepler will actually measure the sum of the planetary and stellar fluxes, 
which means subtracting an estimate of the stellar flux from the observation, and then 
dividing, in order to approach h. 

To complicate matters, h has a zero-point issue, because the minimum of the planetary 
flux — the pedestal of signal on which the variable part of the planetary light curve stands — 
is indistinguishable from steady starlight. Currently, we think the proper way to address 
this issue is to use the zero-mean planetary signal, h-<h>, when comparing theory and 
observation in the process of seeking solutions of the PhO. 

In the ministudy, we assumed that a = 20 ppm applies to measurements of h-<h>, 
and that the photometric errors are normally distributed. 



6. The Merit Function, chiSqR 

For a data set {t, h}, we measure the merit of a candidate solution {a, e, M , u, i, R c s} 
by the sum of the squared, normalized deviations — the "reduced chi square." chiSqR is a 
numerical function with two sets of arguments, the solution and the observational data: 

/(ft [t k ] - <h>) - (h k - <h>) 
chiSqR |{a, e, M , u, i, R eS } , |t, /iJJ = ^ - 1 



k=i 



(20) 



If a global minimum of chiSqR can be discovered, its coordinates are the "solution" — the 
best estimate of the PhO parameters. More generally, surfaces of constant merit in parameter 
space divide more probable solutions from less probable, which is the principle that defines 
confidence regions. If the noise in a Keplerian data set is understood, it may be possible 
to associate absolute levels of probability with such surfaces and regions. In particular, if 
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the measurement errors are normally distributed, and if we assume the Keplerian theory is 
correctly applied, then the probability that the true solution is contained within a surface of 
constant merit can be computed or looked up in a table: it is the cumulative probability of 
the theoretical x 2 statistic. 



7. PhO Fits 

To explore the scientific utility of estimating the six PhO parameters from Kepler data 
sets, we investigated the expected errors in the fitted parameters {a, e, Mq, lo, i, R e s}- We 
used the Levenberg-Marquardt algorithm to locate minima of chiSqR for a subset of the 
Monte Carlo planets in the ministudy. The selection criterion for this subset was pgPeak > 
2.0 x 1CT 8 , which is true for 50% of detected planets and 45% of all planets specified by 
Table 1. 

To save time searching for propitious starting points, the initial estimates of the PhO 
parameters were the true solutions. 

Figures 5-10 show the paramagrams of the parameter-specific median projected normal- 
ized errors (PSMPNEs): 



PSMPNE = Median 



fittedValue — true Value 
maximum Value — minimum Value 



(21) 



where the denominator is the full range of a parameter from Table 1. Table 2 gives the 
median projected normalized errors (MPNEs) and median absolute errors (MAEs) for all 
parameters, for all 10 5 planets considered in this experiment, without regard to location in 
parameter space. 

Figure 8, the paramagrams for the PSMEs of iv, consider the universe of well-inclined 
eccentric orbits (e and i significantly greater than 0). For ~ | or we are looking at 
the planet's fullest brightside or fullest darkside, respectively, at periastron, which amplifies 
light-curve variations that are already most extreme at that time. The distinctiveness of 
this effect at those geometries reduces the errors in u, which is signified by the horizontal 
blue stripes in frames 3, 7 and 10, and vertical blue stripes in frames 13 and 14. Because 
the "backside" light curve changes more rapidly with changes in u, the stripes at u ~ ^ are 
bluer (lower error) than those at uj ~ |. 

Consider the PSMPNEs for i and R e g in Figures 9 and 10. We recognize that performing 
the fits with the observable h-<h> introduces degeneracy in the form of reciprocity between 
Res and i, which will be more problematic for smaller e and upper-mid-range i and for 
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Fig. 5. — The parameter-specific median projected normalized errors (PSMPNEs) for a. Any 
cell in any grid shows the median error in the subject parameter for the subset of a simulated 
data sets with pgPeak > 1.5 x 10 -8 and with true parameters falling in that cell. The color key 
is provided by the histogram,which reports the distribution of projected normalized errors 
(PNEs) by log value. The cyan "J" indicates the body-twin of Jupiter studied in §8. 
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noisier data sets. In this situation, the value of peakToTrough — and approximately the 
whole light curve — can be preserved by increasing i and appropriately decreasing R ef£ , or 
vice versa. This correlation means % and R e ^ are individually less well constrained — i.e., will 
show higher errors — in this regime (see both frames 8). When % is smaller, the effect of its 
variations is larger, and it is better constrained (shows smaller errors); the reverse is true for 
errors in i? e ff, which must fully account for the peak brightness as i — ► | (compare the two 
frames 11). For less noisy data and/or larger e, the information carried in the shape of the 
light curve — i.e., the departures from sinusoidality due to $[/?] and/or e — helps reduce the 
uncertainties due to this degeneracy. 



8. Jupiter Body-Twin 

To examine one PhO in depth, and in order to establish a benchmark planet for the new 
methods and to demonstrate the basic steps for analyzing real data, we used the procedure 
in §3 to create a simulated data set for a Jupiter body-twin, in an eccentric, inclined, 1.9-day 
orbit. The PhO parameters for this data set ("true solution") were 

Ptmc = {a, e, M , u, i, R e s} (22) 
= {0.0300 AU, 0.300, 3.00 rad, 5. 14159 rad, 0.300 rad, 7.90 R®} , 

for which peakToTrough = 2.7 x 10~ 5 . The observational parameters and stellar mass were 
the same as in §3. The data set, with the mean subtracted and folded on the period, is 
shown is Figure 2. 

Figure 4 indicates that the periodogram will detect P truc with 100% certainty. Neverthe- 
less, the signal is less than one-third the median value for detected planets in the population 
of Table 1, and therefore significantly weaker than for any of the planets in the wholesale 
study of §7. Still, this data set constrains the PhO parameters. 

Using the protocol of §4, we find that the periodogram accurately estimates the period 
at T = 1.899 days, which corresponds to a = 0.03001 AU (see Figure 11). 

To find a good starting point to search for a global minimum of chiSqR, we evaluated 
chiSqR directly at a large number of semi-random points in PhO parameter space: the 
semimajor axis was held constant at a = 0.0300 AU, as determined by the periodogram, while 
{M , u, i, R e fi} were drawn from uniform random deviates over their full ranges (Table 1), 
and e was drawn from a uniform random deviate over the range, < e < 0.8. The lowest 
value of chiSqR discovered was 
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{chiSqR direct min , P d irect,min} = {0.992 1 19, {0.03000 AU, 0.300987, 3.04587 rad, (23) 

5.13684 rad, 0.301315 rad, 7.92696 . 

Using indirect, mm as the starting point for the fitting routine, we found 

{ chiSqRfr, min , Pfit.min} = {0.991922, {0.03000 AU, 0.28517, 3.11597 rad, (24) 

5.14031 rad, 0.29315 rad, 8.27574 R e }} , 

which is identical to the fitting result using Ptrue as the starting point. The value of 
chiSqRfn,mm is consistent with a good fit, a valid theory, and correctly estimated, normal er- 
rors. We would expect to encounter a larger value of chiSqR&t^min some 64% of the time. (On 
the same point, the distribution of chiSqR for the fits of §7 closely matched the theoretical 
reduced x 2 probability distribution with 4040 degrees of freedom.) 

-Pfit.min is the best estimate of the PhO. It deviates from the true values by 

P&t mm-Ptruc = {2.9 x 10~ 7 AU, -0.015,0.116 rad, -0.00128 rad, -0.00685 rad, 0.376 , 

(25) 

and the fractional deviation in sini is 2%. For real data, -Ptrue is not known, of course, and we 
are interested in the confidence regions of the parameters, which we expect will be complex 
for any marginal Keplerian data sets (Brown 2004b). That is, we expect degeneracies to 
produce confidence regions with correlations between parameters, and even various types of 
solutions to be compatible with a data set. Such degenerate solutions will evanesce with 
improving signal-to-noise ratio as more information is gathered. Each analysis of a real data 
set is ultimately hand work — teasing out the full information in a data set and finding the 
optimal manner in which to express and apply it. 

Because of the special interest in recovering the inclination angle for computing true 
masses of radial- velocity companions, we computed the confidence regions for i alone, shown 
in Figure 12. Here, we have fixed i to values defined by the abscissa, and achieved the best 
fit to the data by optimizing the other five PhO parameters, then plotted the variation of 
unreduced chi square (chiSq) from the minimum value 

A chiSq = chiSq - 4040 chiSqR^ min . (26) 

AchiSq is distributed as x 2 with 1 degree of freedom, which is the same as the probability 
distribution of the square of a single normally distributed quantity. As shown in Figure 12, 
the confidence region with AchiSq < n 2 can be labeled "n-sigma." In the current case, with 
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Table 2. MPNE and MAE are the median projected normalized 

and absolute errors 



Parameter 


MPNE 


MAE 


Units 


Semimajor axis (a) 


4 x I(T 6 


2 x 1 0~ 7 


AU 


Eccentricity (e) 


0.01 


0.005 




Initial mean anomaly (M ) 


0.003 


0.02 


radians 


Argument of periastron (u) 


0.005 


0.03 


radians 


Inclination angle (i) 
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Fig. 11. — Periodogram of the simulated data set created for the Jupiter body- twin. The 
peak at T = 1.899 days robustly discovers the planet and estimates the semimajor axis 
a = 0.03001 AU with an error of only 0.03%. 
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Fig. 12. — Illustration of sin i recovery for the synthetic data set for a Jupiter body-twin. The 
confidence region for 1-sigma (68%) uncertainty in i is less than ±10%. (Over the indicated 
range, sinz differs from % by less than 3%.) The confidence region becomes asymmetrical at 
higher confidence levels, demanding a thorough treatment of correlated errors. 
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1-sigma confidence (68%), we expect the true value of i to lie in the range 0.265-0.320 rad, 
but with 2-sigma confidence (95%), we can only say it is less than 0.35 rad. (The results 
for sini are basically the same, because % and sini differ by less than 3% over the range of 
Figure 12.) 

9. Phase Function Issue 

The handling of the phase function is an issue. All the results in this paper — except 
in the following paragraphs — used the Lambert phase function in preparing and analyzing 
simulated data sets. Since the phase function of an SPGP is not known in advance, the 
question is, what are the effects of analyzing a data set using the wrong phase function? We 
explore this question by looking closely at the case of the Jupiter body-twin in §8. 

In Figure 3, we introduce a phase-function parameter £ to control a generalized phase 
function ranging continuously from lunar (£ = —1) to Lambertian (£ = +1). This range 
encompasses the variety of phase functions in the Solar System (Sudarsky et al. 2005, Fig. 3). 
At any given phase angle f3, our generalized phase function interpolates linearly between the 
values at phase angle f3 of the lunar and Lambert phase functions. 

It is apparent from Figure 13 that the estimated value of R c g will increase if we use 
for orbit fitting a phase function with a lower value of £ than the true value. Lowering £ 
reduces theoretical values for $, which must be compensated by larger theoretical values of 
R c g. Furthermore, Eq. 1 suggests that the increase in R e g should be proportional to the 
change in the square root of the effective phase integral, /[£], where 

/[£]= f\m(3m,(3]Q[(3]d(3 , (27) 
Jo 

where Q\J3] is the probability distribution of the phase angle (Q\J3] is different for each orbit). 
If we have a good PhO fit using an arbitrarily selected phase function £i, we could roughly 
predict the fitted value of R c s if we had used any other phase function, £ 2 , by multiplying 

R«[ti] by y/IMWS- 

The solid curve in Figure 14 shows the fractional systematic error in -R e ff[£] when the 
Jupiter body-twin data set is analyzed using phase functions in the range — 1 < £ < +1. 
Figure 15 shows Q[/3] for the Jupiter body-twin orbit. The dashed curve in Eq. 14 shows 
the fractional deviation of -R e fr [^] \/ /[^i]/ /[^2] from -R e ff[£i]- This correction indicates a basic 
understanding of the dominant factor in the phase-function issue in regards to R e Q. 

Figure 16 shows the systematic errors due to the phase function for the other five 
orbital parameters in the case of the Jupiter body- twin data set. The semimajor axis a is 
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Phase angle 08) 

Fig. 13. — The generalized phase function developed for investigating the consequences of 
using the wrong phase function when estimating orbital parameters. £ is the phase-function 
parameter. The lunar phase function (£ = —1) is taken from Figure 3 in Sudarsky et al. 
(2005). The Lambert phase function (£ = +1) is given in Eq. 3. For intermediate values of 
£, the phase function at any phase angle f3 is a linear interpolation between the lunar and 
Lambert values at (3. 

, — i — , — , — , — , — , — , — , — , — , — i — , — , — , — , — , — , — , — , — , — i — 




Phase-function parameter (f ) 



Fig. 14. — Solid line: the fractional systematic error in R e g when the Jupiter body-twin 
data set (§8), which was prepared using the Lambert phase function, is analyzed with the 
generalized phase function with parameter £. Dashed line: remaining fractional error after 
the estimate of R e s is corrected by the square root of the change in the effective phase 
integral. 
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Fig. 15. — Empirical probability distribution of the phase angle j3 for the Jupiter body-twin 
data set studied in §8. To obtain this curve, we evaluated (5 at one million random epochs, 
prepared an empirical cumulative probability distribution, took the numerical derivative, 
and verified the normalization to unity of the resulting Q[0\. We used Q[f3] to evaluate the 
effective phase integral for the Jupiter body-twin according to Eq. 27. 
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Phase-function parameter ) 

Fig. 16. — Systematic errors in the PhO parameters a, e, M , ou, and % when the Jupiter 
body-twin data set is analyzed using a phase function different than the Lambert phase 
function (£ = 1) with which the data set was prepared. In the most important case of i, the 
systematic errors are smaller than the typical random errors for similar PhOs (see Fig. 9). 
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not affected, because it is controlled by the periodicity of the signal. The systematic errors 
in e, M , and u are somewhat larger than the typical random errors of similar PhOs (see 
Figs. 6-8. For the inclination angle, i, which is perhaps the most important PhO parameter, 
the systematic error is smaller than the typical random errors: <2% compared with typically 
-10% (see Fig. 9). 

For a real data set, a practical approach would be simply to fit the data using the best 
guess as to phase function, and to qualify the result with that choice. Because the phase 
function itself can only be estimated from data with much higher S/N than afforded by 
Kepler, we will need to be content for now with estimates of e, Mo, £, and R e & that depend 
on the assumed phase function. The most complete treatment will specify an a priori 
probability distribution for the phase function (i.e., for the phase-function parameter £), 
in which case the combined error for each orbital parameter could be expressed as the 
convolution of the probability distributions of systematic and random errors. 

Because of the somewhat smaller systematic error in inclination angle i as compared 
to its random errors, our prediction that Kepler-level photometry may break the msini 
degeneracy in a useful manner is not substantially undermined by the phase-function issue. 

10. Discussion and Conclusion 

We have defined and put forward an approach to investigating the PhO, which we 
believe will come to be recognized as a Keplerian entity on a par with radial-velocity and 
astrometric orbits. 

The PhO solves in principle one of the outstanding problems in astronomy: the m sin % 
degeneracy in radial-velocity observations. The way is now open to computing true masses 
for radial- velocity companions, starting with SPGPs. 

Brown (2009) shows for large, noisy Keplerian data sets — particularly in the case of 
astrometry, but radial velocity and photometry are all basically the same in this regard — 
that search completeness via periodogram, as well as the completeness of accurate estimates 
of orbital parameters by least squares, are solely functions of the signal-to-noise ratio (S/N). 
S/N is defined as the semi-amplitude of the signal times the square root of the number of 
data points, divided by the single-measurement error. For example, 50% search completeness 
calls for S/N = 5-8, due to any combination of signal amplitude, number of data points, and 
single-measurement error (so long as the orbit is well covered by the measurements). (The 
color scale above the abscissa in Figure 4 shows S/N fa 7 for the PhO.) High completeness 
for, say, ±10% accuracy in estimates of orbital parameters typically calls for S/N = 50-60. 
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In the case of the Jupiter body-twin data set shown in Figure 2 and treated in §8, S/N = 48. 
Therefore, detection is guaranteed for this planet for the assumed Kepler performance and 
observing parameters. Furthermore, the high S/N of the Jupiter body-twin data set means 
that the small deviations of the fitted orbital parameters from the true values (Eq. 25) are 
no surprise. 

PhO observations may be difficult. There is uncertainty about how low the albedos of 
SPGPs actually are, about their phase functions, and about a possible noise floor due to 
systematic errors from instrumental and stellar sources. Nevertheless, simple detection of 
SPGPs in reflected light should be robust in the regime of Kepler photometry, and estimates 
of all six orbital parameters may be feasible in at least a subset of cases. 

Burrows et al. (2008) offers the most recent discussion of observed and theoretical values 
of geometric albedo p for the case of SPGPs. The best measured value, upper limit p = 
0.038 ± 0.045 of Rowe et al. (2008), from MOST observations of HD 209458b, is compatible 
with cloudless models. The same models suggest an that p could increase by factor two 
at the somewhat longer wavelengths sampled by Kepler. Meanwhile, transit photometry of 
SPGPs has revealed astonishingly large radii — 1.74 -Rj up i t er in the case of TrES-4 (Mandushev 
et al. 2007). A planet with the size of TrES-4 and twice the geometric albedo of HD 209458b 
would offer R c s = 5.4 ± 5.8i? ea rth- If that planet were on the orbit of the Jupiter body-twin, 
the signal-to-noise ratio would be S/N = 22+^- 

As expected for any new method, particularly one debuted with such epochal claims as 
breaking the msini degeneracy, our treatment of the PhO will certainly attract questions 
and concerns, which we must study and address. We would be the first to caution about 
systematic effects, particularly about a noise floor that might limit the averaging out of 
random measurement errors. Nevertheless, the potential payoff of determining the inclination 
angles of radial-velocity companions independently — that alone — is so great that we must 
press forward. Whatever its current technical limitations, we can expect the grasp of this 
new method to improve with time, as instrumentation improves and systematic errors are 
addressed. 

A key issue is how well our method stands up to more realistic noise, particularly due to 
stellar and planetary variability. The latter must wait for Kepler results and the possibility 
that observations may discover variations in planetary properties. The former, however, is a 
known issue, which we must study using more realistic random deviates for the photometry. 
That research will likely produce selection criteria for the future target lists. We expect 
that a meaningful subset of Kepler targets will support PhO modeling. This expectation is 
supported by the prediction of Jenkins & Doyle (2003) that Kepler will discover "from 100 
to ~760 SPGPs" in reflected light by periodogram searching. All these SPGPs are potential 
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candidates for PhO analysis. 

It has been asked how our simple model and six-parameter solutions relate to the ap- 
parent complexities of the forward problem, that is, predicting photometric properties from 
detailed planetary models (Seager et al. 2000; Dyudina et al. 2005; Sudarsky et al. 2005; 
Gaidos et al. 2006). We are observers, working on the reverse problem, which is inferring 
planetary properties from a photometric record. Unless the signal-to-noise ratio of Kepler 
data sets is much higher than we expect, we predict satisfactory x 2 probabilities for fits 
using only six parameters, which would say that additional parameters are not required by 
the observations. For systematic effects with assumed values that are not constrained by 
the observations, the associated range of uncertainty must be convolved with random errors 
when the total uncertainty in estimates are expressed, as we have outlined in §9 the case of 
the phase function. 

We thank Adam Burrows for his comments on the manuscript. We thank Christian 
Lallo for aiding the computations and assisting with the graphics. We thank Sharon Toolan 
for her expertise preparing the manuscript. JPL contract 1289452 with the Space Telescope 
Science Institute provided support for this research. 
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